Energy density functional study of nuclear matrix elements for neutrinoless (3(3 decay 
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We present an extensive study of nuclear matrix elements (NME) for the neutrinoless double 
beta decay of the nuclei 48 Ca, 76 Ge, 82 Se, 96 Zr, 100 Mo, 116 Cd, 124 Sn, 128 Te, 130 Te, 136 Xe, and 150 Nd 
based on state-of-the-art energy density functional methods using the Gogny D1S functional. Beyond 
mean-field effects are included within the generating coordinate method with particle number and 
angular momentum projection for both initial and final ground states. We obtain a rather constant 
value for the NME's around 4.7 with the exception of 48 Ca and 150 Nd, where smaller values are 
found. We analyze the role of deformation and pairing in the evaluation of the NME and present 
detailed results for the decay of 150 Nd. 

PACS numbers: 21.60.Jz. 23.40.-s, 23.40.Hc 



Double beta decay is an extremely rare process where 
an even-even nucleus decays into the even-even neighbor 
bypassing the energetically forbidden odd-odd interme- 
diate isobar. Along the nuclear chart, there are only few 
candidates, all found in the valley of stability [1]. The 
double beta decay where two electrons and two neutri- 
nos are emitted in the final state {2v(5(3) has been ob- 
served experimentally in several isotopes with half-lives 
~ 10 19-21 yr. This process conserves the leptonic num- 
ber and is compatible with Majorana or Dirac neutrinos. 
There is also a second mode without neutrino emission 
(CV/3/3) that is possible only if the neutrinos are massive 
Majorana particles and is related to the absolute mass 
scale of these elementary particles [1]. Except to one 
controversial claim [2] , 0^/3/3 decay has not been detected 
and is currently the main goal of several projects world- 
wide [3] . The half-life of this process between + states 
for mother and granddaughter nuclei can be written as 

[1]: 
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where (mpp) is the effective Majorana neutrino mass, 
m e is the electron mass, Goi is a kinematical phase 
space factor, and finally, M 0l/ is the nuclear matrix el- 
ement (NME). Eq. 1 shows that a precise determination 
of the effective neutrino mass requires, apart from the 
experimental measurement of the 0vj3j3 half-life, a re- 
liable calculation of the NME. So far, several nuclear 
structure methods have been employed. The most used 
among them are the Interacting Shell Model (ISM) [4, 
5], proton- neutron Quasi- Random Phase Approximation 
(QRPA) [6-9], and, more recently, Projected Hartree- 
Fock-Bogoliubov (PHFB) in limited configuration spaces 
and using schematic Pairing plus Quadrupolc interac- 
tions [10, 11] and Interacting Boson Model (IBM) [12]. 
In this letter, we present major improvements with re- 



spect to previous PHFB calculations [10, 11]: we use 
state-of-the-art density functional methods based on the 
well established Gogny D1S functional [13] and a much 
larger single particle basis (eleven major oscillator shells); 
we perform particle number and angular momentum pro- 
jection for both mother and granddaughter nuclei; and 
include configuration mixing within the generating co- 
ordinate method (GCM) framework [14, 15]. All these 
developments have been shown to be necessary for an 
unified description of nuclear structure [14] . In particu- 
lar, particle number projection before variation is funda- 
mental for a correct treatment of pairing correlations [15] 
that play a very important role in 0v/3/3 decay [4]. 

The Qv(3(3 NME can be separated into three terms: 
Fermi (F), Gamow-Teller (GT) and Tensor (T) [1]: 
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with gv = 1 and gA = 1.25. The tensor term is small 
according to the ISM and QRPA calculations [5, 8] and 
will be neglected in this work. We use the closure ap- 
proximation [1] to sum over intermediate states in the 
odd-odd nucleus as currently it is not possible to com- 
pute odd-odd nuclei using beyond mean-field methods: 
symmetry restoration methods including blocking effects 
have not been fully developed so far. Therefore, a cal- 
culation of the 2v(3j3 mode cannot be performed as the 
closure approximation is not valid in this process. The 
different terms in Eq. 2 can be expressed as the expec- 
tation value of a two-body operator between the initial 
and final states, i.e. M°!j GT = (0^|M°^|0+). Detailed 

expressions for M^f,Q T can be found in Ref . [5] . We have 
included high order currents [6], nucleon finite size cor- 
rections [6] and radial short range correlations treated 
within the Unitary Correlator Method [8, 16]. In the 
GCM+PNAMP approach (shortly GCM from now on), 
the initial (i) and final (/) many-body wave functions 
are found as linear combinations of particle number N, Z 
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and angular momentum 1 = projected wave functions 
with different intrinsic quadrupole deformations /3: 

\0 + )=J29pP I= °P N P Z \^) (3) 
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where P I=0 , P , P z are the corresponding angular mo- 
mentum (/ = 0) and particle number projectors [17]. 
The intrinsic axial symmetric Hartree-Fock-Bogoliubov 
(HFB) wave functions |$^) are solutions to the vari- 
ation after particle number projection equations con- 
strained to a given value of the axial quadrupole defor- 
mation, /? [15, 18]. Therefore, intrinsic deformation of 
the system is naturally included in the formalism and 
pairing correlations properly taken into account. Fi- 



nally, the coefficients g@ are found by solving the Hill- 
Wheeler- Griffin (HWG) equation [17, 19]. First, for each 
nucleus we transform the non-orthogonal set of wave 
functions {P I=0 P N P z \$p) } into an orthonormal one 

{l A > = E/3 ^W= pI=0pNpZ \®p)} b y diagonalizing the 
norm overlap matrix, J2(i'(^p\P I=0 P N P Z \^P') u A,p' = 
n-AUA,i3- In this basis, the HWG equation reads: 
J2a> £ AA'G A , = E a G%, where e\A' are the so-called en- 
ergy kernel [15, 19]. Finally, the coefficients for the lowest 
eigenvalue are used to compute both the so-called col- 
lective wave functions F((3) = J2a Ga u a,/3 - probability 
distribution for the state to have a given deformation - 
and the NME: 
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Particle number conservation together with using large 
and identical configuration spaces for protons and neu- 
trons guarantees that Ikeda's sum rule is fulfilled. Addi- 
tionally, to our knowledge, this is the first implementa- 
tion of the GCM method for calculating transitions be- 
tween different nuclei including particle number symme- 
try restoration in both the initial and final states. 

We now present our results for the 0^/3/3 NME dis- 
cussing in detail the decay of 150 Nd. To check the re- 
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FIG. 1: (color online) Comparison between theoretical and 
experimental ground state bands for 150 Nd and 150 Sm. Inset: 
collective wave functions of the ground states as a function 
of the intrinsic deformation. The values for 150 Sm have been 
shifted up by 0.05. 



liability of the method for describing properties of the 
initial and final nuclei, we show in Fig. 1 a comparison 
between the experimental and theoretical ground state 
bands for 150 Nd and 150 Sm. We observe a rather good 
agreement for both excitation energies and £?(E2) transi- 
tion probabilities, with the theoretical results predicting 
a slightly larger rotational (collective) character than the 
experiment. The computed double beta decay Q- value 
is 2.99 MeV while the experimental value is 3.37 MeV. 
The inset of Fig. 1 shows the probability distribution 
for the mother and granddaughter states to have a 
given intrinsic quadrupole deformation j3. Both 150 Nd 
and 150 Sm have well deformed prolate ground states with 
P ~ +0.40 and j3 sa +0.25, respectively. These values 
are compatible with the rotational bands shown in the 
figure. Equation 4, shows that the Qv(3j3 NME can be 
expressed as a sum of matrix elements between states 
of different intrinsic quadrupole deformation for the ini- 
tial and final nuclei. Consequently, the matrix element 

<*„ \P N I P Z I M^P'^P^P^^^) 
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plicitly the strength of the Qvf3j3 operators as a func- 
tion of the deformation of the nuclei involved in the de- 
cay. In Fig. 2(a) we show these matrix elements for GT 
component (the Fermi part gives a similar pattern and 
it is not shown). The strength is concentrated rather 
symmetrically in the diagonal part of the figure, imply- 
ing that the decay between states with different initial 
and final deformation is hindered. Moreover, spherical 
configurations are the most preferred to decay. Never- 
theless, Fig. 2(a) shows also an interval of deformation 
close to the spherical configuration where non-negligible 
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FIG. 2: (Color online) (a) Strength of the GT operator as a 
function of the intrinsic deformation of the mother and grand- 
daughter nuclei. Shaded area corresponds to the maximum 
probability for finding the mother and granddaughter states, 
(b) Pairing energies as a function of the intrinsic deformation. 



off-diagonal matrix elements (greater that 0.5) are ob- 
tained (f3 G [—0.2, +0.2]). The absolute maximum is 
found at (+0.03,0) with additional local maxima, for 
example at (+0.52, +0.52). These results are in agree- 
ment with ISM and PHFB calculations that have shown 
a significant decrease of the NME with increasing differ- 
ence in quadrupole deformation between initial and final 
states [10, 20]. This trend has also been observed in 2v(3(3 
QRPA calculations [21]. The final value of the NME is 
determined by weighting the strength of the transition 
operator with the wave functions of the initial and final 
states which selects the relevant region of deformation. 
This area is marked by a shaded circle in Fig. 2(a) show- 
ing that in this case both the difference in deformation 
and the absence of shape mixing inhibit the transition. 
The final result for the NME is M°» = 1.71 with 1.28 
and 0.43 coming from Gamow- Teller and Fermi compo- 
nents, respectively. In order to shed additional light on 
the structure of the strength we show in Fig. 2(b) the 
pairing energy —E pp of the nuclei involved in the decay. 
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FIG. 3: (Color online) Nuclear matrix elements calculated 
for different methods (ISM [5, 22], QRPA(Jy) [8], QRPA(Tu) 
[7], IBM-2 [12], PHFB [10]) with UCOM short range correla- 
tions. QRPA values are calculated with gA = 1-25 and IBM-2 
and PHFB results are multiplied by 1.18 to account for the 
difference between Jastrow and UCOM [23]. 



This energy is defined as the pairing tensor part [14, 17] 
within the PNAMP approach. We find a strong correla- 
tion between the structure of the NME and —E pp , both 
having maxima at similar values of the deformation for 
mother and granddaughter nuclei. This result is also in 
agreement with ISM and QRPA calculations where the 
largest values for the NME are obtained with well paired 
wave functions with large zero seniority components in 
the spherical basis [4, 7]. 

We now present the values of the NME (Eq. 4) com- 
puted for several double-beta emitters and compare with 
the ones calculated with other methods that use similar 
assumptions concerning the neutrino potentials. In Fig. 
3 we observe that the NME's are rather constant around 
an averaged value of 4.7 for the decay of 76 Ge, 82 Se, 96 Zr, 

100 MOj n6 Cdj 124 Sn! 128 Te; t30 Te and tSG^ bemg 96 Zr 

the one with the largest value. In these nuclei, the differ- 
ences in deformation between the initial and final states 
are not very significant and also the structure and ab- 
solute values of the transition strength are quite similar. 
There are two exceptions to this general trend, namely, 
48 Ca and 150 Nd where the NME is significantly smaller. 
In the later, the difference between deformation of the 
mother and granddaughter is remarkable (inset Fig. 1) 
and this is precisely the main source of suppression of the 
NME. On the other hand, due to the double magic char- 
acter of 48 Ca, this nucleus has the smallest value for the 
pairing energy and transition strength among the nuclei 
studied in this work. Concerning the NME's obtained 
by other methods, we attribute first the small NMEs ob- 
tained by PHFB to the reduced configuration space used 
and the lack of particle number restoration and shape 
mixing in those calculations. We also show that the low- 
est NMEs are obtained within the ISM, a rather con- 
stant value around 2.5 except for 48 Ca. We expect that 
enlarging the model space used in the ISM calculations 
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will enhance the values of their NMEs [4, 9, 22]. Con- 
sidering higher seniority components in our method by 
including quasiparticle excitations in the intrinsic wave 
functions may reduce slightly our values. Some differ- 
ences are also found between the GCM and QRPA val- 
ues. The main differences between these approaches are 
the assumption of spherical symmetry in QRPA, the ab- 
sence of quasiparticle excitations in the GCM approach 
and the much larger -and no core- single particle basis 
used in GCM. Furthermore, neither triaxial, mirror, nor 
time reversal symmetry breaking effects are included in 
our GCM calculations because they are beyond current 
computational capabilities. We expect that the inclusion 
of these degrees of freedom will not change significantly 
the structure of the ground states as the nuclei studied 
here are either spherical or well deformed. To validate 
our approach we have computed the total GT strengths 
£+(-) for the decay of granddaughter (mother) nuclei 
defined, for example, in Ref. [30]. In Table I, we com- 
pare the calculated £+(-) with the experimental values 
measured in charge exchange reactions. As in QRPA 
and ISM calculations, a quenching factor of (0.74) 2 has 
been introduced [21, 30]. Finally, we evaluate the half-life 
of each nuclei based on the NME calculated with GCM 
method. In Table I, we show the difference between the 
calculated and experimental Q-values and we observe an 
excellent agreement in most of the cases except for 96 Zr 
and 100 Mo, where an over binding of 96 Mo and 100 Ru iso- 
topes gives such differences. These are precisely the de- 
cays with largest NME and smallest half-lives predicted 
by our calculations. 
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In summary, we have presented the first calculations of 
Oz^/3/3 decay within the energy density functional frame- 
work including beyond mean-field effects. We have an- 
alyzed the role of the intrinsic quadrupole deformation 
and pairing content of the nuclei involved in this pro- 
cess. Decays between spherical initial and final shapes 
are found to be favored while large differences in deforma- 
tion hinder significantly the transition probability. Our 
calculations constitute the first consistent evaluation of 
the Oi//3/3 decay of 150 Nd. 
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TABLE I: Difference between theoretical and experimental 
Q values, total GT strength for granddaughter (mother) S+ 
(S-) , NME and predicted half-lives for several decaying 
nuclei assuming {mpp) = 0.5 eV. 
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